# -*- coding: utf-8 -*-

# -- Sheet --

import numpy as np


def Gauss_CPE(A, b):
    n = len(b);
    index = 1;
    x = np.zeros(n)
    for k in range(n):
        a_max = 0
        for i in range(k, n):
            if abs(A[i][k]) > a_max:
                a_max = abs(A[i][k])
                r = i
        if a_max < 0.00000001:
            index = 0
            return
        if r > k:
            for j in range(k, n):
                z = A[k][j];
                A[k][j] = A[r][j];
                A[r][j] = z
            z = b[k];
            b[k] = b[r];
            b[r] = z
        for i in range(k + 1, n):
            m = A[i][k] / A[k][k]
            for j in range(k + 1, n):
                A[i][j] = A[i][j] - m * A[k][j]
            b[i] = b[i] - m * b[k]
    if abs(A[n - 1][n - 1]) < 0.0000001:
        index = 0
        return
    for k in range(n - 1, 0 - 1, -1):
        for j in range(k + 1, n):
            b[k] = b[k] - A[k][j] * x[j]
        x[k] = b[k] / A[k][k]
    return index, x


def Jacobi(A, b, x0, it_max, ep):
    D = np.diag(np.diag(A))
    U = -np.triu(A, 1);
    L = -np.tril(A, -1)
    B = np.dot(np.linalg.inv(D), (L + U))
    f = np.dot(np.linalg.inv(D), b)
    x = np.dot(B, x0) + f
    k = 1;
    index = 1
    while np.linalg.norm(x - x0) >= ep:
        x0 = x
        x = np.dot(B, x0) + f
        k = k + 1
        if k > it_max:
            index = 0;
            break
    return k, index, x


def G_S(A, b, x0, it_max, ep):
    D = np.diag(np.diag(A))
    U = -np.triu(A, 1);
    L = -np.tril(A, -1)
    B = np.dot(np.linalg.inv(D - L), U)
    f = np.dot(np.linalg.inv(D - L), b)
    x = np.dot(B, x0) + f
    k = 1;
    index = 1
    while np.linalg.norm(x - x0) >= ep:
        x0 = x
        x = np.dot(B, x0) + f
        k = k + 1
        if k > it_max:
            index = 0;
            break
    return k, index, x


A = [[15, 3, 5, 6], [2, 15, 6, 5], [3, 6, 18, 4], [2, 5, 3, 13]]
b = [22, 28, 13, 15]
x0 = [0, 0, 0, 0]

print(Gauss_CPE(A, b))
print(Jacobi(A, b, x0, 100, 0.000001))
print(G_S(A, b, x0, 100, 0.000001))
